#This to obtain:
#Figure_3b.jpg

library("foreign")
library("ggplot2")
library("sf")
library("RColorBrewer") # creates nice color schemes
library("scales")       # customize scales
library("gridExtra")    # mutiple plots
library("plyr")     
library("dplyr")
library("grid")
library("tidyr")
library("ggsn")
library("lemon")
library("reshape2")
library("lfe")
library("broom")
library("readxl")


#setwd("C:/Users/bogdanp/")
setwd("/Users/bgpopescu/")

#Reading Polygons
HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")

HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))
HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)


final_sp<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))
data<-subset(final_sp, bosnia_border_mun==0)

m1 <- felm(pct_no_water~treat+
             POINT_X+POINT_Y+zagreb_distance+
             bfe1+bfe2+bfe3+bfe4+
             bfe5+bfe6+bfe7+
             quad1+ quad2 
           |0|0|NAME_1, 
           data = data)
summary(m1)
m1<-tidy(m1)

m2 <- felm(pct_no_sewer~treat+
             POINT_X+POINT_Y+zagreb_distance+
             bfe1+bfe2+bfe3+bfe4+
             bfe5+bfe6+bfe7+
             quad1+ quad2
           |0|0|NAME_1, 
           data = data)
summary(m2)
m2<-tidy(m2)

m3 <- felm(pov_rate_income~treat+
             POINT_X+POINT_Y+zagreb_distance+
             bfe1+bfe2+bfe3+bfe4+
             bfe5+bfe6+bfe7+
             quad1+ quad2 
           |0|0|NAME_1,
           data = data)
summary(m3)
m3<-tidy(m3)

#This saves all the previous regression results in a dataframe.
all_models <- bind_rows(
  m1 %>% mutate(model = "pct_no_water"),
  m2 %>% mutate(model = "pct_no_sewer"),
  m3 %>% mutate(model = "pov_rate_income"))

#This saves the OLS mode by gathering variables.
ols_table <- all_models %>%
  dplyr::select(-std.error, -statistic, -p.value) %>%
  gather(key, value, estimate)

#The variables become column names.
#Models become rows.
regcoef<-dcast(ols_table, model~term, fill=0)
colnames(regcoef)
colnames(regcoef)[2] <- "Intercept"

########
#Part 2#
########

#### Load data that will be used in all the figures

boundary <- st_read(dsn='./Dropbox/Legacies_Project/Analysis/data/data.gdb',
                   layer="krajna_line_HRV")
grid <- st_read(dsn='./Dropbox/Legacies_Project/Analysis/data/data.gdb',
                   layer="grid4c")
boundary <- st_transform(boundary, crs = st_crs(grid))

grid_df<-subset(grid, select = c("xcoord", "ycoord", "treat", "NEAR_FID", "dist_krajna6", "dist_zagreb"))
grid_df <- grid_df%>% mutate(zagreb_dist = (dist_zagreb/1000),
                             krajina6_dist = (dist_krajna6/1000))

grid_df$quad1<-ifelse((grid_df$xcoord > 18 & grid_df$xcoord<20) | 
         (grid_df$xcoord<47 & grid_df$xcoord>44), 1, 0)

grid_df$quad2<-ifelse((grid_df$xcoord > 16 & grid_df$xcoord<18) | 
                        (grid_df$xcoord<47 & grid_df$xcoord>44), 1, 0)

grid_df$quad3<-ifelse((grid_df$xcoord > 14 & grid_df$xcoord<16) |
                        (grid_df$xcoord<47 & grid_df$xcoord>44), 1, 0)


#The next line is to get border FE
grid_df$krajina6_id1<-0
grid_df$krajina6_id1[grid_df$krajina6_id==1]<-1
grid_df$krajina6_id2<-0
grid_df$krajina6_id2[grid_df$krajina6_id==2]<-1
grid_df$krajina6_id3<-0
grid_df$krajina6_id3[grid_df$krajina6_id==3]<-1
grid_df$krajina6_id4<-0
grid_df$krajina6_id4[grid_df$krajina6_id==4]<-1
grid_df$krajina6_id5<-0
grid_df$krajina6_id5[grid_df$krajina6_id==5]<-1
grid_df$krajina6_id6<-0
grid_df$krajina6_id6[grid_df$krajina6_id==6]<-1
grid_df$krajina6_id7<-0
grid_df$krajina6_id7[grid_df$krajina6_id==7]<-1
grid_df$krajina6_id8<-0
grid_df$krajina6_id8[grid_df$krajina6_id==8]<-1

regcoef[is.na(regcoef)] <- 0


##############
#pct_no_water#
##############

subgrid <- grid_df
subdata <- data
subdata <- subset(subdata, krajna6_distance<=30)
subcoef <- subset(regcoef, model=="pct_no_water")
subgrid <- (subgrid  
            %>% mutate(ones=1))  # Create a column of ones for constant
subgrid <- data.frame(lapply(subgrid, function(x) as.numeric(as.character(x))))



subgrid$fit1<-subcoef$treat * subgrid$treat+
  subcoef$POINT_X * subgrid$xcoord+
  subcoef$POINT_Y * subgrid$ycoord+
  subcoef$zagreb_distance * subgrid$zagreb_dist+
  subcoef$bfe1 * subgrid$krajina6_id1+
  subcoef$bfe2 * subgrid$krajina6_id2+
  subcoef$bfe3 * subgrid$krajina6_id3+
  subcoef$bfe4 * subgrid$krajina6_id4+
  subcoef$bfe5 * subgrid$krajina6_id5+
  subcoef$bfe6 * subgrid$krajina6_id6+
  subcoef$bfe7 * subgrid$krajina6_id7+
  subcoef$quad1 * subgrid$quad1+
  subcoef$quad2 * subgrid$quad2+
  subcoef$Intercept * subgrid$ones



pct_no_water<-ggplot()+
  geom_raster(data=subgrid, aes(x =xcoord, y = ycoord, fill=fit1))+
  geom_sf(data = boundary, colour = "red", linewidth = 1.5)+
  geom_point(data=subdata, aes_string(x='POINT_X', y='POINT_Y', fill="pct_no_water"), size=1.25, 
             colour="gray50", shape=21, na.rm=TRUE)+ #pch=21 allows for borders
  scale_x_continuous("Longitude") +
  scale_y_continuous("Latitude") +
  theme_bw()+
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  scale_fill_continuous(high="black",low="white", name="Percent")+
  ggsn::scalebar(x.min = 18.5, x.max = 19.1,
           #Above you mention how long the scalebar should be
           y.min = 46.2, y.max = 46.4,
           #Above you mention how thick the scalebar should be
           dist = 30, dist_unit = "km",
           transform = T, model = "WGS84",
           location = "topright",
           st.size = 4,
           #Above you have the font size of the numbers below the scalebar
           st.dist =0.4,
           #Above you have the distance between the bar and the text, as a proportion of the y axis.
           height=0.2)
pct_no_water<-reposition_legend(pct_no_water, 'bottom right')
pct_no_water

ggsave(pct_no_water,file="./Dropbox/Legacies_Project/Paper/figures/Figure_3b.jpg", 
       height=10, width=20, 
       units = "cm", dpi=300)

